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NUMERICAL COMPUTATION OF EXPONENTIAL MATRICES 
USING THE CAYLEY-HAMILTON THEOREM 


Haivey Walden 

NASA/Goddard Space Flight Center 
E. C. Roelof 

The Johns Hopkins University 
Applied Physics Laboratory 

ABSTRACT 

A method for computing exponential matrices, which often arise 
naturally in the solution of systems r ^ linear differential equations, is 
developed. An exponential matrix is generated as a linear combination of a 
finite number (equal to the matrix order) of matrices, the coefficients of 
which are scalar infinite sums. The method can be generalized t o apply to 
any formal power series of matrices. In this paper, attention is focused 
upon the exponential function, and the matrix exponent is assumed tri- 
diagonal in form. In such case, the terms in the coefficient infinite sums 
can be extracted, as recursion relations, from the characteristic polynomial 
of the matrix exponent. Two numerical examples are presented in some 
detail: (1) the three-dimensional infinitesimal rotation rate matrix, which is 
skew-symmetric, and (2) an N-dimensional tri-diagonal and symmetric 
finite difference matrix which arises in the numerical solution of the heat 
conduction partial differential equation. In the second example, the known 
eigenvalues and eigenvectors of the finite difference matrix permit an analyt- 
ical solution for the exponential matrix, through the theory of diagonaliza- 
tion and similarity transformations, which is used for independent verification. 
The convergence properties of the scalar infinite summations are investigated 


V 
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for finite difference matrices of various orders up to ten, and it is found that 
the number of terms required for convergence increases slowly with the order 
of the matrix. Successive squarings of the exponential matrix demonstrate 
that the resultant limiting matrix tends toward zero in a smooth fashion, due 
to the dominance of the largest eigenvalue of the finite difference matrix. 

An appendix to the paper presents an algorithm for the evaluation of expo- 
nential matrices. 
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NUMERICAL COMPUTATION OF EXPONENTIAL MATRICES 
USING THE CAYLEY-HAMILTON THEOREM 


INTRODUCTION 

The concept of exponential matrix arises naturally in physical problems in which the under- 
lying (possibly high-order) differential equations can be reduced to a system of linear differential 
equations. In the simplest case, the fundamental solution y(x) of the linear differential equation 
dy/dx = Ay, where A is a scalar constant, is given by: 

y(x) = y(0) e^^ = y(0) S A^ x*^ , 

k=o k! 


When y is a vector and A is a square matrix, then the infinite series of matrices given above is still 
well defined, i.e.^, converges uniformly and absolutely for x in any bounded interval [ 1 , pp. 64-65] , 
and provides the solution to the linear differential equation in terms of the exponential matrix, 
exp (Ax). 

As an example of a physical problem in which exponential matrices arise, consider the repre- 
sentation of the propagation operator for energetic charged particles in the interplanetary magnetic 
field. The propagation equation, a second-order partial differential equation with three independ- 
ent variables, may be reduced, under appropriate assumptions, to a steady-state relation with only 
two independent variables. This steady-state second-order partial differential equation may be 
solved by finite-difference methods using discrete ordinates, in which the second-order partial dif- 
ferential operator is replaced by a second-order difference operator represented by a matrix M, 
which is tri-diagonal in form. The propagation equation may thus be written as 


aWj 

as' 


S M. W. 
j=i J 
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whose solution may be expressed in terms of an exponential matrix; 

W(s') = W(0) e^s' = ^(0) 2 (s')^ . 

k»0 k! 

A method for generating such an exponential matrix as a linear combination of a finite num- 
ber N of matrices will be discussed in this paper. 

POWER SERIES OF MATRICES 

According to the Cayley-Hamilton theorem [1, p. 113] , a square matrix satisfies its own 
characteristic equation. For an N x N matrix M, let the characteristic polynomial be written as 

p'^(X) = det(M„ -M) = I X" , ( 1 ) 

where the coefficient aj^ = 1 , and Ij^ is the N x N identity matrix. Thus, by the Cayley-Hamilton 
theorem, p^ (M) = 0, or 

MN = _ ^ M" . (2) 

This relationship leads to the conclusion that any podtive integral power of M can be expressed as 
a linear combination of the N matrices, Ij^ , M, , . . ., . In fact, when the inverse M~^ exists, 

the preceding statement may be generalized to any integral power. Suppose this is true for some 
integral power k, so that 

k N-l 1, n 

M*" = - A|; M" , (3) 

where the a|^ are suitably chosen scalar coefficients. Then, multiplying by M, one finds 

T.,k+1 Ak xk //IN 

M = - 2^ Aj^_j M -- Ajq_j M . (4) 

Substituting for from equation (2), 
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M 


k+1 _ N-l 


\k X . ,n j. *k „N , 

n?i ^^N-1 ^ ^N-1 > 


( 5 ) 


whereas equation (3) for k+1 yields 


M' 


k+1 _ Ak+1 wH A k+1 1 

- \ M - Ap In 


( 6 ) 


Equating coefficients of equal powers of M in equations (5) and (6) results in the recursion relation- 
ships: 


= ^n -1 “ An_i aJJ, forl<n<N-l, 


.k+1 ^ .k N 
Aq ” An _i a^ . 


(7a) 

(7b) 


k _ 


If we establish the convention A_j = 0, then the above relations combine to 


A 


k+1 _ Ak 


k _N 


An-1 - An_i for 0<n<N-l . 


( 8 ) 


By comparing equation (2) with equation (3) for k - N, we find 


An ~ ^n ’ 0<n<N-l 


(9) 


In what follows, we shall assume that k > N in equation (8) and extend the range of the indices of 
a]^ by defining (in terms of the Kronecker delta) 


An==-\n» forO<k<N-l. 


( 10 ) 


Any formal power series f(M) derived from the infinite expansion, 


f<X) = S x'‘ , 

k=0 ^ 


( 11 ) 


can be expressed as a linear combination of the N matrices, In , M, M‘‘ , . . , in the form 
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f(M) = S Bk 

k=0 ^ 


S 
k = 0 


N-1 
Bu S 
n»0 




(12a) 


N-1 

s 

n*=0 


oo 

( S 
k=0 




(12b) 


upon substitution of equation (3) and interchanging the order of the summations. Note that eval- 
uation of f(M) in this manner requires only N scalar infinite sums as coefficients of the matrices. 

A typical power series might be the exponential function f(X) = e^ for which Bj^ = 1/k! . In 
this case, the sequence Bj^ is rapidly converging toward zero, and, in order to overcome possible 
numerical difficulties in a situation in which the coefficients Aj^ are diverging as k-^oo, we may 
define a new well-bounded sequence 


C 


k 

n 


= 



(13) 


For this power series, equation (12b) may be written 


M N-1 

n = 0 


OO 

( S 

k=0 




(14) 


in which the left side is an exponential matrix. The recursion relation (8), upon multiplication by 
\+\> may be rewritten as 

“ “^N-l “n ). for 0<n<N-l,andk>N, (15) 

- 1 k 

since Bj^^ j /Bj^ = (k+l)“ . The convention C_j = 0 is required, as well as the analogues of the 
initialization equations (9) and (10), which are 


, N 

C;; = forO<n<N-l, (16) 

y ’ forO<k<N-l. (17) 
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Note that this method requires explicit computation of tlie inverse factorials only up to the matrix 
order N, in contrast to the usage of equation (12b). 

The remaining task is the extraction of the coefficients aj^ from the characteristic polynomial, 
det (XIj^ - M), For a general NxN matrix, there are coefficients my involved, and the task is 
algebraically complicated. However, for a tri-diagonal matrix M, the coefficients aj^ may be deter- 
mined recursively. First note from the determinant 


P (X) = 


X-iHj, 

-mi2 

0 

0 • • 

• • 

. . . . 

. . . . 0 

-m2i 

X—m22 


0 . • 

• • 

. . , . 

. . . . 0 

0 

• 

• 

• 

-m3 2 

• 

• 

• 

• 

• 

• 

• 

-m34 ■ 

• 

• 

• 

• 

• 

• • • • 

. . • • 0 

■ 

• 

• 

• 

• 

0 . 

f . . • • 

• 

• 

• 

• 

• 

■2 

• 

• 

• 

1,N-1 

• 

• • 
• 

0 . 


. . . . 0 


-mw 

,N-1 

^-•nN.N 


(18) 


the relationship obtained by expanding up the last column; 


P^ (X) = (^) -1 ,N "’n ,n -1 P^ > 


N-2, 


(19) 


where p^**^ and ^ are the appropriate sub-determinants of p^ (X). Substituting the definition 
(1) in equation (19), we obtain 


N 


N ^n 


2 a-X“ = 

n=0 “ 


V .n 

a_ , a 


n=l 


‘n-1 


m 




N.N “n 


m» 


N-1,N ”^N,N-1 “n 




( 20 ) 


Equating coefficients of like powers of X on both sides, we obtain the recursion relations: 


“ ^N-i “ ^ definition) , 


(21a) 
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N-l 

“M-2 


- JTljj , 


(21b) 


a 


N 

n 



N-l N-2 

%,N ~ "’n-I.N "^N,N-1 » 

for 1 <n <N— 2, 


(21c) 


-m 

”^N,N ^0 


*N-1,N 


N-2 

*N.N-1 “O 


(21d) 


The recursion relations (21) clearly also hold for any NxN submatrix of M containing the elements 
my where 1 < i, j < N, so that directly we have: 


p\X) = A-m 


11 


P^(X) 


-mu 


“«^21 

X-ltlj2 


= 'hr ~ (mjj +m22)X + (m^j mj2 - nij, ni2i). 


(22a) 


(22b) 


Since p^(X) = aj + aj X and p^(X) = aj f aj X + a^ X^, it is seen that the initial values for 
the recursions are 

aj = 1; aj = -m^^ , (23a) 

^2 ~ ~~ (^11 "^ ^22^ ’ ^0 ~ ^11 ^22 ~ ^12 ^21 * (23b) 

The coefficients a^j may be generated recursively by using the formulas (2 1) with N replaced by N 
for N = 3 to N = N. Alternatively, equation (19) will generate p (X) from p (X) and p (X) = 1, so 
the recursive technique can actually start with N = 2, 


6 
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numerical example 1; ROTATION RATE MATRIX 

As an example, consider the 3x3 infinitesimal rotation rate matrix [2, p. 127] , which happens 
to be skew-symmetric: 


CO = 


0 cOj — CO 2 

-CO 3 0 co^ 

CO 2 — coj 0 


(24) 


If we set CO 2 = 0, then co is a tri-diagonal matrix, and the techniques of the previous section may be 
applied, The initial values (23) are 


4 = -"’ll = 0- 

aj - "(1^11 + ^22^ ^ 

^0 ~ ^22 ~ ^12 ^21 “ ^3 ’ 

and the recursion relations (21) yield 

al = -m33 al - aj = ° » 

= ^0 - '^33 - ^23 ^32 a{ = CO 3 + COj = \u?\^ , 

4 = - *^33 = 0 • 

Thus, equation (2) becomes 


(25a) 

(25b) 


(25c) 


(26a) 

(26b) 

(26c) 


co' = - I a'co" 

n=0 " 


= —jCOf CO 


(27) 


At this point, the tri-diagonal restriction (viz., that co, = 0) can be relaxed, since the result (27) also 
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2 2 2 

holds for coj ^ 0, with the gent'ializc-d definition that |w| s + coj + <*>3 , as may be verified by 
direct calculation. The recursion iclations (8) become 

= Ajj_i -Ajlcol^ej^, for 0<n<2 and k>3, (28) 

since aj = Icol^ 8 j for 0 < ii < 2 as a result of equations (26). Evaluation of equation (28) 


produces the following results: 

K = <29a) 

K = K-l - At l“l' «i„ = -l“l" »i„ . (29b) 

aJ = A'„.i - A* lup = -Icol" S2„ , (29c) 

a’ = a‘_i -A*|co|^Si„ = (29d) 

A* = a’„.i - a51co1'8,„ = |co 1®82„. (29e) 
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The function 


f(cjt) = S 
k «0 


( 32 ) 


becomes, upon substitution of equation (30), 


f(«t) = B,!, + w 


k =0 


. r V n / i\k 1 i2k .2k+2, 2 

+ [ S (“1) loj| t ] CO 

k =0 ^ 


For the function f(X) = e^ for which = 1/k! , this becomes 


(33) 


e“' = I 3 + 


r “ (~icoi^ t ~i . Ff (-itop t^)*‘ tH 

Lic=o ( 2 k+l)! J “ [k=o ( 2 k+ 2 )! J 


CO 


(34) 


or, upon recognition of the bracketed terms as MacLaurin trigonometric series expansions for the 
sine and cosine, 


cot , . / sin lco| t \ 

" “ '3 "TT" " ^ 


/ 1 - COS |coj t \ 

V / 


(O^ 


(35) 


Note that if Icoj > 1 , then the calculation of the a|J coefficients, as shown in equations (29), pro- 
duces an unbounded sequence. This may readily be avoided by the use of the recursion (15) and 

Ic Ic lc 

the calculation of the coefficients of definition (13) in place of where the sequence of Cj^ 

converges to zero. Specific numerical values of these coefficients are shown in Table 1, for the 

rotation rate matrix with co^ = 1, CO 2 = 0 , and CO 3 - 3. 

In the calculation of the infinite summations appearing in equation (14), the relative error in the 
partial sums after each block of N terms is computed and compared to a pre-selected convergence 
criterion e in order to evaluate numerically whether convergence has been attained. That is, con- 
vergence is deemed to occur when 


9 
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PN-I u (p-l)N-l ,, 

s q - s 

k=0 " k=0 ” 


pN-1 

S 

k = 0 



< e for 0 < n < N-1 and p = 2, 3, , . . (36) 


Partial sums are compared only after each block of N terms is accumulated into the summation in 

order to rule out false indications of convergence due to a given term cjj vanishing. For a rather 

stringent convergence criterion of e = 1 0“^ the summation S c!f converged after 33 terms 

k=o ^ 

(p= 11 blocks) to the value ~sinlcj|/lcj| s 6.5407 x 10"^, and the summation S C« converged 

k=o ^ 

after 30 terms (p= 10 blocks) to the value -(1 -cos |co|)/|(op = -0.19998, where |cop = 10 
for the matrix previously cited. 


Table 1 

Comparison of Selected Values for the Scalar Coefficients 
(Rotation Rate Matrix with = I, = 0, cOg = 3) 


Indices 

Scalar Coefficients 

n 

k 

aJ; [eq. (8)] 

cj; [eq.(15)J 

1 

1 

-1 

-1 

1 

5 

-100 

-0.833 

1 

15 

10^ 

7.65 X 10-^ 

1 

25 

-10^2 

-6.45 X 10-1^ 

1 

35 

1012 

9.68 X 10"24 

1 

45 

O 

7 

-8.36 X 10“25 

2 

2 

-1 

-0.5 

2 

6 

-100 

-0.139 

2 

10 

-10^ 

-2.76 X 10-2 

2 

20 

10^ 

4.11 X 10-10 

2 

30 

-lOi-l 

-3.77 X 10-19 

2 

40 

1019 

1.23 X 10-29 
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NUMERICAL EXAMPLE 2: FINITE DIFFERENCE MATRIX 

As a second example, consider the NxN tri-diagonal and symmetric matrix: 


T = 


-2 1 


1 -2 


0 0 • • • 0 


1 0 • • • 0 


0 1 -2 1 • • • 0 


0 0 0 0 • • • 1 


0 0 0 0 • • • 0 


0 0 


0 0 


0 0 


1 -2 


where the diagonal entries are all —2 and the super- and sub-diagonal entries are all unity. This 
matrix arises frequently in the numerical solution of partial differential equations, in particular in 
the explicit finite difference approximation to the heat conduction equation [3, pp. 60-64} . The 
eigenvalues of Tj^ are 




and the corresponding eigenvectors are 


rnr 2nir Nnir 

sin . sm , . . . , sin — — 

N+1 ’ N+1 ’ ' N+1 


0 ‘ 


for n = 1, 2, . . N, 


where the superscript “t” indicates the transpose. These values can be verified by substitution into 


T* V. u^ — * Um . 

N n n n 


The initial values (23) are, in this case, 


1 ^ ^ A 2 o 


and the recursion relations (21) yield 
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.N OaN-1 _ a^-2 

^0 ^0 


N ^ N-i ^ N-1 
®n ^n-l ^ “An 


aN = a^“^ + 2 , 
Sn-I “N-2 ’ 


aIJ"^ forl<n<N-2, 


(41a) 

(41b) 

(41c) 


so that, with N replaced by N “ 3, we find 


®o " " ^0 " 


aj = ag + 2aJ - a[ = 10 


= a^ + 2 — 6 , 


(42a) 

(42b) 

(42c) 


and, with N replaced by N = 4, we find 


aj = ag + 2aJ — aj = 20 , 

^2 " ^^2 " ^2 ” » 
a3 = aj + 2 = 8 . 


(43a) 

(43b) 

(43c) 

(43d) 


N 

This recursion process may be continued, and Table 2 displays the full array of coefficients aj^ for 

matrix order N = 8, The final row of Table 2 thus provides the coefficients in the characteristic 

8 2 
polynomial p (X) of equation (1), as well as the coefficients of each of the powers Ig , Ig , Tg , . . 

Tg appearing in the summation of equation (2). Upon calculation, each of the sequences for 

n = 0, 1 , 2, . . ., 7 of equation (1 5) is seen to converge toward zero as k -► «». In the calculation of 

the infinite summations appearing in equation (14), convergence is based on the inequality (36) for 
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Table 2 


Full Array of Coefficients aJJ 
(Finite Difference Matrix of Order N = S) 



the same criterion of e = 10"^ ^ as used in the previous example. Results of applying this conver- 
gence criterion are shown in Table 3. It is noted that convergence in all 8 summations is attained 


Table 3 

oo 1^ 

Convergence Properties of the Summations S C„ 

k=o " 

(Finite Difference Matrix of Order N = 8 with Criterion e = 10”^ 


Summation 

n 

Convergence after 

Approximate 
Summation Value 

k terms 

p blocks 

0 

40 

5 

-0.99996 

1 

40 

5 

-0.99944 

2 

40 

5 

-0.49779 

3 

40 

5 

-0.16272 

4 

48 

6 

-3.7853 X 10-2 

5 

48 

6 

-6.1599 x 10-2 

6 

48 

6 

-6.2691 X 10-^ 

7 

48 

6 

-2.9581 X 10-^ 
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after either 5 or 6 blocks of 8 terms each and that the converged absolute values of the summations 
monotonically decrease with the power of Tg with which they are associated in equation (14). 

The quantity exp (Tj^ ) may be computed numerically by equation (14). With the known 
eigenvalues (38) and eigenvectors (39) of Tj^ , the exponential matrix may readily be determined 
analytically as well. Since the N eigenvalues of Tj^ are distinct, Tj^ is similar to a diagonal matrix 
A = diag (Xj , X 2 , . . ., Xj^ ) whose entries are the eigenvalues of Tj^ [ 1 , p. 77 1 , or 

Tj, = SAS"‘ , (44) 


where the similarity transformation matrix S is formed of ordered columns composed of the eigen- 
vectors of Tj^ . Furthermore, since Tj^ is symmetric, its eigenvectors are mutually orthogonal [ 1, 
p. 99] . If, in addition, the eigenvectors (39) are normalized, so that Uj^' ^ ^ ® simi- 

larity matrix S' formed of ordered columns composed of the Uj^' vectors will be a unitary matrix 
such that (S')~^ = (S')^ By the “telescoping” property of the similarity transformation, for any 
positive integer k. 


= (SAS‘^)^ = (SAS'^)(SAS"^).,.(SAS”^) fork factors 

= sa’^s"^ , 


(45) 


so that 


exp (Tj^ ) 


00 

S 

k = 0 



«* SA^S"^ 

S — 77 ^" 

k=o k! 



- S [exp(A)j S'^ . 


Since A is a diagonal matrix, A^ - diag (X^ , X^ > • • . » x|^), and exp (A) = diag (e^^ , e^ 
If, in equation (46), S' is substituted for S, then the final result is 


(46) 
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expdj^) = S^ldiag(e 


\2 




)! (S')‘ . 


(47) 


For the value N = 8 , Table 4 presents the eigenvalues (38) and unnormalized eigenvectors (39) of 
Tjj . where, due to the periodic properties of the sine function, there are only four non-zero magni- 
tudes involved in the 64 eigenvector components. Note that the eigenvalues span more than an 


Table 4 


Eigenvalues and Unnormalized Eigenvectors 
(Finite Difference Matrix of Order N = 8 ) 


Eigenvalues 
\ [Eq. (38)] 



Eigenvector Components^ 

[Eq. (39)1 



«2 

^3 

^4 

U 5 


u? 

«8 


= -0.12061476 



^3 

S4 

S4 

S3 

S2 

Si 


= -0.46791111 

^2 

S4 

S3 

Si 

-Si 

-S3 

-S4 

-S2 

^3 



S3 

0 

-S3 

-S3 

0 

S3 

S3 


= -1.65270365 

S4 

Si 

-S3 

-S2 

S2 

S3 

-Si 

-S4 


= -2.34729636 


-Si 

-S3 

S2 

h 

-S3 

-Si 

S4 


= -3 

S3 

-S3 

0 

S3 

-S3 

0 

S3 

-S3 

X? 

= -3.53208889 


-S4 

S3 

-Si 

-Si 

S3 

-S4 

S2 

^8 

= -3.87938524 

Si 

-§2 

S3 

-S4 

S4 

-S3 

S2 

-Si 


♦where, by definition, Sj = 



sin 7 r /9 = 0.34202014, 
sin27r/9 = 0.64278761, 
sin rr/3 = VJ/2 = 0.86602540, 
sin 47 r /9 = 0.98480775. 


order of magnitude, since the maximum ratio is \/\ — 32. Since the eigenvectors are all of mag- 
2 2 

nitude |u| = |Uj^| == 9/2 for n = 1, 2, , . . , 8 , the normalized similarity transformation matrix S' 

is given explicitly by 
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where = sin tt/9, S2 = sin 27 t/9, S3 = sin tt/ 3 =>/J/2, and s^ = sin 47 t/ 9. Note that S' is sym- 
metric as well as unitary, so that (S')^ = (S')“^ = S'. On inserting equation (48) into equation 
(47), it is found that exp (Tg) is both symmetrical and “back-symmetrical”; that is, the exponen- 
tial matrix exhibits symmetry about both principal diagonals. Explicitly, the exponential matrix is 
of the form 


exp(Tg) = 




®3 

®4 



®7 


^2 


®10 

^11 

®12 

®13 

«14 



®10 

«15 

®16 

®17 

®18 

®13 




®16 

®19 

®20 

®17 

®12 



®12 

®17 

®20 

®19 

®16 


®4 


®13 

®18 

®17 

®16 

^15 

®10 

®3 


®J,4 

®13 

®12 


^‘'10 

®9 

®2 



^6 

65 

®4 

®3 

®2 



(49) 


where both the analytical and numerical values for each of the 20 distinct elements ej are given in 
Table 5. 
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Table 5 

Analytical and Numerical Values for Exponential Matrix Elements 
(Finite Difference Matrix of Order N = 8) 


Symbol 

Analytical Expression* 

— 
Numerical Value 


+ 44 + 4e^6 + 44 

0.21526929 

®2 

^ 1 ^ 2^18 ^ 2 ^ 4^27 ^ S ^36 ® 1 ® 4^45 

0.18647807 


® 1 ® 3^18 ^ 2 ^ 3^27 ® 3 ^ 4^45 

8.6.3736679 x 10"^ 

«4 

^ 1 ^ 4^18 ® 1 ^ 2^27 " ® 3^36 ® 2 ® 4^45 

2.74614615X 10"^ 


® 1 ® 4^18 ~ ® 1 ^ 2^27 ^ ® 3^36 “^ 2^4 ^45 

6.64880517 X 10"^ 

®6 

® 1 ® 3^18 “ ^ 2^®27 ^ 3 ^ 4^45 

1.29935583 X 10‘^ 

®7 

® 1 ^ 2^18 “ * 2 ^ 4^27 S ^36 “ ^ 1 ^ 4^45 

2.12770688 X 10""^ 


® 1^18 ~ ^ 2^27 ^ 3^36 ® 4^45 

2.95813145 X 10"^ 


2 1 2 ¥- 1 *^ 1 _2 1 _2 TQ *^ 

^ 2^18 ^ 4^27 ® 3^36 ® 1®45 

0.30164296 

®10 

^ 2 ^ 3^18 ® 3 ® 4®'27 ® 1 ^ 3^45 

0.21393953 

®11 

^ 2 ^ 4^18 '*' ® 1 ® 4^27 "■ ^ 3^36 “ ^ 1 ^ 2^45 

9.30224731 x 10“^ 

®12 

^ 2 ^ 4^18 ~ ® 1 ® 4^27 ” ® 3^36 ^ 1 ^ 2^45 

2.87608174 X 10"^ 

^3 

^ 2 ® 3^18 ■“ ^ 3 ^ 4^27 ^ 1 ^ 3^45 

6.86157586 X 10“^ 

®14 

® 2^18 ~ ^ 4^27 ® 3®36 "“ ®1 ^45 

1.32893715 x 10"^ 

®15 

4<4 + 4 + 4 ) 

0.30829176 

®16 

® 3 ^ 4^18 ■‘' ^ lS ^27 ^ 2 ^ 3^45 

0.21523888 

®17 

^ 3 ^ 4^18 “ ^ lS ®'27 " ^ 2 ^ 3^45 

9.32352438 x 10"^ 

®18 

4®T8 - e;, - E«) 

2.87903987 x 10”^ 

^19 

^ 4^18 ’*' ^ 1^27 ® 3®36 ® 2^45 

0.30850453 

®20 

* 4^18 ~ ^ 1^27 S ^36 ■" ® 2®45 

0.21526847 


•where, by definition, Ej^g=e^^± e^®, E^^=e^^ ± , 

Egg =gX4 + gXs 
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OO 

The convergence properties of the infinite summations S appearing in equation (14) 

k«o " 

were further investigated for finite difference matrices of orders N = 2, 4, 6, 8, and 10. The results 
are summarized in Table 6 for the convergence criterion e - 10”^^ used earlier as well as’ for a more 
relaxed criterion of e ~ 1 0”^^ . It is seen that use of the more stringent criterion leads to a modest 


Table 6 


OO 

Convergence Properties* of the Summations S 

k=o 


for Finite Difference Matrix of Order N 


A. Using Criterion e = 1 


Summation, n 
Matrix Order, N 

Number of terms, k required for convergence 

n 

N 

0 

1 

2 

3 

H 

5 

6 

B 

8 

9 

2 

24 

24 









4 

28 

28 

28 

32 







6 

30 


36 

36 

36 

36 





8 

32 

32 

40 

40 

40 

40 


40 



10 

40 

40 

40 

40 

40 

40 

40 

40 

40 

40 



*The number of blocks, p required for convergence is given by p = k/N. 
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increase in tlie number of terms required for convergence, Note also that the number, p = k/N of 
blocks required for convergence actually decreases as the matrix order N increases: for N = 2, p = 

12 or 15; for N ~ 4, p 7, 8, or 9; for N = 6, p ~ 5, 6, or 7; for N = 8, p = 4, 5, or 6; and for N = 

10, p - 4 or 5. The summations denoted by the lower integer values n converge to larger absolute 
values than those denoted by larger values of n (see Table 3), and the former converge more rapidly, 
viz., generally, convergence is attained one block earlier. The most significant result indicated by 
Table 6, however, is that as the order N of the matrix increases, the number k of terms required for 
convergence of the infinite summations increases very slowly. 

A final investigation was made into the properties of the exponential matrix raised to an inte- 

j? 

gral power. Specifically , calculations were made to evaluate the matrices (exp Tj^ ) , where the 
exponent J2 is itself a positive integral power of two. Table 7 presents the range of matrix elements 
for the matrices (exp Tg )^, where 12 = 2^, for q = 0, 1 , . . . , 9, These matrices are readily evaluated 


Table 7 


Range* of Elements for Successive Squarings of Exponential Matrix of Order N = 8 


Value of 
Power, C 

Minimum Element of 
Matrix (expTg)® 

Maximum Element of 
Matrix (exp Tg)^ 

1 

2.96 X 10"^ 

3.09 X 10-^ 

0 

6.82 X 10“'^ 

2.07 X 10-^ 

4 

4.70 X 10-3 

1.40 X 10-^ 

8 

7.79 X 10-3 

8.28 X 10-3 

16 

3.72 X 10-3 

3.13 X 10-3 

32 

5.48 X 10-^ 

4,54 X 10-3 

64 

1.15 X 10-3 

9.57 X 10-3 

128 

5.13 X 10-9 

4.25 X 10-8 

256 

1.01 X 10-^3 

8.39 X 10-^3 

512 

3,94 X 10-39 

3.26 X 10-38 


*A11 elements of all matrices are positive. 


19 










by successive squarings of the form (exp Tg ^ -• (exp Tg ^ ^ (exp Tg ^ \ so that for £ = 2^, 
only q successive matrix squarings are necessary. An analytical check is readily at hand, sinct by 
the "telescoping” property of the similarity transfonriation, both sides of equation (47) c'.n be 
raised to the power £ as 

[expert)]® = S'[diag(e®’^‘,e'’^2,...,e'’'N)] (S')*, (50) 

Table 7 shows that as £ increases, the maximum elenjent of (exp Tg monotonically decreases, 

The minimum element of (exp Tg)^ initially Increases with £, but then later also decreases as £ 
increases. Most significantly, the range of elements within the matrix (exp Tg )® monotonically 
decreases as £ increases and (exp Tg approaches the zero matrix as £ -►‘» in a rather smooth 
fashion. The smoothness of the convergence toward zero, and the decreases in the range of the 
matrix elements, are due to the dominance of the largest eigenvalue, Xj (refer to Table 4). 
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APPENDIX 

ALGORITHM POR EVALUATION OF EXPONENTIAL MATRICES 

1 . The following parameters are given initially: 

a. the matrix order N, where N > 2 is a positive integer; 

b. the upper limit, on the number of recursions calculated, where > N is a 

positive integer; 

c. the criterion, e, for convergence of the partial sums, where 0 < e < I ; and 

d. the tri-diagonal square matrix M of order N, This can be efficiently provided by ini- 
tializing M as an N X N zero matrix and then substituting the elemsents on the principal 
diagonal, the super-diagonal, and the sub-diagonal, as follows: 

mil for i = 1,2,. . „N, 

*^i,i+l . . .,N-1, 

mi+i i for i = 1,2,.. ., N-1. 

N 

2. The initial values for the array are: 

4 = -"'ll 

- "'ll *"22 — *"12 *"21 (23)* 

aj = -mjj - mj2 

a" = 1 for n = 0, 1,..,,N (21a) 

3. The array values a^ are calculated recursively for N - 3, 4, . . ., N from: 

N _ N-1 ^ S-2 

^0 ””^N,N % " *^N-1,N "^N,N-1 % 

a” = aJJ,! - - ”^n-i,n “’n.n-i ^ii n= 1, 2, . , . ,N~2 

N _ N-1 _ 

^N-1 - “n_2 " "^N,N 

*These equation numbers provide a reference to the discussion in the main bcxiy of the text. 
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4. TJie initial values for the array Cjj are; 
N 

d = ^ for n « 0, 

" N! 


c|5 = 0 if k n 
c!; = -1/k! if k = n 


font, k - 0, 1, . . . , N-1 


5, The array values C]J are calculated recursively for k = N, N+1, . . . , from; 


'^ N-1 ‘*0 


'0 

ik+l 

'n 


k+1 


‘t (cJ;., - cjj., a|J) for n = l,2...„N-l 


k+1 


06) 

(17) 


(15) 


6. Evaluate integral powers of the tri-diagonal matrix M iteratively, as follows: 
M® = In . the NxN identity matrix, 

^ M, 

and M for k = 2, 3, . . . , N~l, 

or, in terms of matrix elements: 

Ic N k-l 
mil = 11 m.f, nu. 

ij c=i ifi fij 


oo ^ 

1. Calculate successive partial sum approximations to the infinite sums £ C„ for n = 0, 1, . . ., 

k=o " 

N— 1, and rompiire the relative differences in successive partial summations after each block of 
N terms to the convergence criterion, where 


(p+DN-l u 
£ 

k=0 


n 


pN-1 I. 

- £ C 

kto " 


(p+DN-l J. 

s < 

k = 0 " 


< e for n = 0, 1, . . . , N~1 


(36) 


indicates convergence of the summation for some value, p = 0, 1 , . . . (By convention, define 
pN-l t, 

£ C„ = 0 for p = 0,) 
k=o " 
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If (be above inequality is not satislied Tor all /em and positive inlegi'a! values ol' p such that 
(pt‘l)N 1 ^ oomeiitence of (be Inilnlle sum has not been attained and either 

^max iUL'ivased. 


8. I'be exjHvnentlal matrix is evaluated as a linear eombmation oi‘ N matrices» as Tollows; 




i-M « ( s o;; ) M 

11 »U k“0 


a 


(b 


Note: The above abiorilUm Tor the evaluation of exponential matriees has been implemeiUed as 
a FORTKAN language computer pix>gram and is available, upon tvquest, IVom the authors. 



